Seismic source identification of the 9 November 2022 Mw 5.5 offshore Adriatic sea (Italy) earthquake from GNSS data and aftershock relocation

The fast individuation and modeling of faults responsible for large earthquakes are fundamental for understanding the evolution of potentially destructive seismic sequences. This is even more challenging in case of buried thrusts located in offshore areas, like those hosting the 9 November 2022 Ml 5.7 (Mw 5.5) and ML 5.2 earthquakes that nucleated along the Apennines compressional front, offshore the northern Adriatic Sea. Available on- and offshore (from hydrocarbon platforms) geodetic observations and seismological data provide robust constraints on the rupture of a 15 km long, ca. 24° SSW-dipping fault patch, consistent with seismic reflection data. Stress increase along unruptured portion of the activated thrust front suggests the potential activation of longer portions of the thrust with higher magnitude earthquake and larger surface faulting. This unpleasant scenario needs to be further investigated, also considering their tsunamigenic potential and possible impact on onshore and offshore human communities and infrastructures.

www.nature.com/scientificreports/ significant coseismic offsets in relative time-series, which sum the cumulative deformation of the two M > 5 seismic events (in Fig. S1 of Auxiliary materials we report the SENI and FANO time series example). Taking advantage of some onshore (storage centers) and offshore (seabed-anchored hydrocarbon platforms) industrial infrastructures hosting continuous GNSS stations, integrated with public and private GNSS networks onshore, we collected a precious dataset to infer additional constraints on the activated buried thrust. This is the first time that GNSS observations collected on offshore infrastructures are used to constrain the coseismic deformation pattern of an earthquake. Offshore platforms generally stand on legs of steel or concrete that are piled into the sea floor, and thanks also to their own weight remain fixed in place, securely anchored to the sea floor. Their stability has been also proven by the long-term GNSS time-series acquired by on-board continuous stations 9,10 therefore they are suitable for active tectonic studies 1 . Here, by processing the extensive dataset collected at both onshore and offshore GNSS stations, we estimated the cumulative coseismic deformation pattern related to the 9 November 2022 seismic events. Despite its millimeter-scale, the coseismic displacement field clearly depicts a reverse faulting mechanism pattern, providing key constraints for the seismic source identification (Fig. 1) and subsequent stress transfer analysis.
In this study, we show a fast preliminary space-time characterization of the seismic sequence obtained by seismic data from the INGV monitoring system and the coseismic slip as constrained by GNSS data, coupled with available subsurface public seismic reflection data. Our main aim was to constrain the source of the 2022 Adriatic earthquake for its use in future studies of different kinds (e.g., hazard mitigation, earthquake and tsunami scenarios, hydrocarbon platform safety).

Figure 1.
Epicentral area of the 9 November 2022 offshore earthquake (yellow star). Black dots are the aftershocks. Green arrows are the horizontal coseismic displacements recorded by the GNSS stations (triangles); tringle colors indicate the vertical coseismic displacement as indicated in legend. We also report the main thrusts 6 close to the seismic sequence and the moment tensor calculated for the mainshock (http:// terre moti. ingv. it/ event/ 33301 831).

Geological setting
In the northern Adriatic Sea, the external portion of the Apennine orogenic belt involves the sedimentary succession of the Adria passive margin and the overlying foredeep deposits associated with the Neogene to recent W-directed Apennines subduction 11 . In this sector of the Apennines orogen, the thrust front is buried beneath Early Pliocene-Quaternary synorogenic deposits, as indicated by seismic-reflection data 6,12-15 . The thrust system is characterized by ramps climbing off the basal detachment, located within the Triassic evaporites, to join upward a shallower detachment located along Early Miocene marls 1,15,16 . Several seismic sources were associated with these structures, as also suggested by instrumental and historical seismicity 5,17 . Recently released GNSS data from offshore hydrocarbon seabed-anchored platforms, located few tens of kilometers from the epicentral area, permitted to define a shortening rate of about 1.5 mm/yr 2,9 mostly accommodated along the Apennines frontal thrusts in a SW-NE direction 1 , significantly improving our knowledge of the shortening rate accommodated between the Northern Apennines and the Adria plate 18,19 .
The frontal sector of the Northern Apennine orogen is characterized by NW-SE striking and NE-verging double-plunging anticlines developing above SW-dipping ramps (Fig. 2). The November 2022 seismic sequence is located along the most external structure of the Apennine orogen [20][21][22] . In the epicentral area, only poor-quality vintage public seismic profiles are available, and it hazardous, to discriminate fault planes at depth (Fig. 2b). We tentatively plotted the time-converted hypocenter along the SV-167-4-83 seismic profile (http:// www. videpi. com) using two end-member constant velocity scenarios (i.e., 4000 m/s and 5000 m/s) to highlight the possible depth range of the mainshock hypocenter (Fig. 2b). Seismic line B-402, located 10 km far from the mainshock and available from ViDEPI public database (http:// www. videpi. com), provides a general characterization of the structural setting of this frontal sector of the Northern Apennine (Fig. 2c). Line B-402 allowed the recognition of a flat-ramp-flat geometry for the external thrust, with the lower flat not clearly detectable but interpreted to be deeper than 4 s TWT (Fig. 2c).Considering the hypocentral depth of the mainshock and aftershocks, we suggest that the activated portion of the thrust corresponds to the ramp, intermediate between the deeper and shallower detachments, and involves the Mesozoic carbonate sequence, likely down to the Triassic evaporites.

Earthquake locations and seismological parameters. A total of 594 aftershocks in the period 22
November 2022-05 January 2023 and the main shock have been relocated by using the hypoDD relative method 8 on P-and S-wave arrival times picked by the seismologists on duty in the INGV control room. Results differ from those of the INGV control room (available at http:// terre moti. ingv. it/): in particular, the M w 5.5 is located about 4.4 km south and 4.1 km deep with respect to the control room one. We estimate the uncertainties on the location parameters through a bootstrap approach (see the Methods section for the details on the location procedure and error estimation). Aftershock relocations (Fig. 3) show a NNW-trending, 25° SSW-dipping, 15 km long fault, seismically defined between 5 and 10 km depth. The original data from the INGV seismic control room have been analyzed to characterize the magnitude distribution and time decay of the aftershocks. We estimated magnitude distribution of the seismic sequence following the Gutenberg-Richter 23 recurrence relation. The inferred magnitude of completeness (Mc) value is rather high (Mc = 2.0), as expected for the large distance to the coast (15-20 km) and the high noise of the closest stations installed on thick sediments of the Adriatic foreland, and the b-value is 0.94 (Fig. 3).
We plot the number of aftershocks/day and the daily average energy release as a function of time, above the completeness magnitude threshold of 2.0, which is reached 4 h after the mainshock (Fig. 3). We applied the corrected Akaike criterion to verify whether the decay rate is better fitted by a regular Omori-like (single slope) model, or by a diffusion-Omori (double slope) model 24 , finding that the two-slope model is statistically better than the single-slope behavior, with a 55% likelihood. In the first six days after the mainshock, the slope of the line (t −0.5 ) has the characteristic decay rate of a diffusion-like phenomenon 25 , whereas after the sixth day the slope agrees with a regular Omori-like decay rate (slope ≈ 1). These observations argue for a substantial contribution of fluid pressure diffusion to the seismogenesis in the early stage (6 days) of the sequence.
Ground displacement and seismic source model. Measured GNSS coseismic displacements well captured the cumulative coseismic deformation of the two M w > 5 earthquakes, with average 1σ uncertainty of ~ 1 and 3.6 mm for the horizontal and vertical components, respectively. The largest horizontal displacements have been observed at the stations located along the coast and shows a centripetal pattern pointing toward the epicentral area (Fig. 1). The 4 closest offshore stations constrain the displacement field in the north and east zone. DARB, the closest station, and FAUZ record about 4 mm northward and about 2 mm westward, respectively. BREN and CLAW show submillimeter displacements, constraining the source length along the strike.
To define the seismic source responsible for the mainshock, we modeled the coseismic displacements using the Okada 26 formalism. Geometry and kinematics of the fault are obtained in an initial unconstrained non-linear inversion of GNSS data. Best solution well reproduces the displacement pattern (Fig. S2), with a WNW-ESE striking and 24° dipping reverse fault plane, consistent with both moment tensor solution (striking 128°) and the aftershocks distribution described in the previous paragraph. Thus, we fixed to this value the strike angle of the geodetic solution. Uncertainties and trade-off parameters are shown in Fig. S3. Coseismic slip distribution along the identified fault plane is calculated by means of linear inversion. Slip along our solution concentrates between 4.5 and 10 km depth, with a maximum of about 12 cm (Fig. 4). Synthetic displacements satisfactorily fit the observed ones at the onshore stations and at on offshore ones located within a radius of 30 km from the epicenter (see also Table S1). Although noisy, the model partially reproduces the observed vector at FAUZ station. Farther ones show no displacements or values close to the uncertainty, constraining the fault dimension. www.nature.com/scientificreports/  www.nature.com/scientificreports/ We also investigated the relationships between the mainshock source and the closest faults by means of the Coulomb Failure Function (hereinafter ΔCFF). We calculated the ΔCFF along the closest faults from DISS catalogue 5 , offshore and along the coast, caused by the 9 November earthquake. Thus, as input, we used the slip distribution estimated during the linear inversion. Stress transfer analysis reveals stress variation lower than the minimum stress-variation-threshold magnitude of 0.01 Mpa 27 along the faults located westward of the seismic source (Fig. S4). Faults along the strike (NW and SE) of the source show values up to 0.07 MPa (to the North), and up to 0.32 Mpa along the unslipped fault portion during the 9 November earthquake.

Discussion
The DD relocation shows that aftershocks align on a ca. 25°, SW-dipping plane between 5 and 10 km depth, although the absolute hypocentral depth of the mainshock suffers from large uncertainty. Our modeling shows a maximum slip of 12 cm on a low-angle thrust fault, whose geometry (strike and dip angles) is consistent with that indicated by aftershocks distribution, confirming the reliability of our modeling results. The coseismic slip mainly concentrates at depth from the mainshock hypocenter, suggesting a downward and coast-ward directivity of the coseismic rupture.
The sequence has a magnitude distribution of events with a b-value close to 1 and shows a first diffusive pattern (0.5 slope in Fig. 2d) followed by a normal decay (slope 1), with an oscillating release of energy (Fig. 2c) that tends to be Omori-like. The aftershocks diffusion is usually interpreted as a diffusion of the stress induced by the mainshock, either by fluid transfer in the crust 25,28 , or due to a viscous relaxation process 29 . Indeed, the seismicity diffusion may also result from a simple cascade process where the mainshock triggers aftershocks that in turn trigger their own aftershocks, leading thus to a general expansion of the aftershock volume. We hypothesize that after an external forcing that locally increased the seismicity rate during the very first part of the sequence, diffusion was reduced leading to a standard decay. This inferred Omori trend is like that of the Emilia 2012 seismic sequence 30 , but the initial diffusive pattern was limited to the first 4-5 days and did not lead to the occurrence of further mainshocks. Trying to estimate in semi real time such patterns could be relevant www.nature.com/scientificreports/ to make inference on the possible evolution of the sequence, but caution should be exercised when using data of an ongoing seismic sequence 31 . The compressional structure (fault-related fold) is well illuminated in public seismic reflection profiles but the definition of the activated thrust portion is hindered by the high noise versus signal ratio affecting the seismic data at depths larger than 2 s TWT (Fig. 2b). Shallow buried anticlines are well mapped in the upper 5-6 crustal kilometers [20][21][22] . They developed above 20°-35° dipping ramps splaying from a regional basal décollement located at depth and dipping between 1° and 7° toward the west 32 . One of the above-mentioned ramps seems to be the fault plane activated during the earthquake and it might correspond to the NAF-10 structure 21 . Similarly to the 20 May 2012 Mw6.0 Emilia and the 26 November 2019 Mw6.4 Durres compressional earthquakes, also in this case the coseismic slip propagated from the hypocenter downward 30,33,34 , a feature typical of thrust faults 35 . This suggests that a crucial role in the future evolution of the sequence could be played by the stress increase (up to 0.31 Mpa) along the un-slipped portions of the activated thrust, and, more in general, along the thrust-front. Under this assumption, we expect to observe an along-strike progression of the sequence with minor impacts on the coastal belt.

Conclusions
The analysis of geodetic and seismological data available since a few days after the 9 November 2022 earthquake in the northern Adriatic area defines a set of constraints about the seismic sequence and criticalities in the source determination, as well as the crucial role of geophysical/geodetic data acquired by hydrocarbon industry infrastructures. Monitoring of buried faults in offshore areas can take great advantage from the presence of industrial infrastructures that could host GNSS continuous stations. This aspect, coupled with amphibious seismic networks are indispensable for monitoring such faults and constrain seismological parameters.
Hypocenters distribution and GNSS coseismic displacements inversion reveal that slip occurred along a 15 km long fault patch along a ca. 24° SSW-dipping thrust fault, consistent with seismic reflection data propagating downward from the mainshock hypocenter. Slip amount rapidly decreased upward and along strike, where stress accumulated. If, on one hand, the downward pattern of slip intuitively reduces the tsunamigenic www.nature.com/scientificreports/ potential, on the other hand the stress increase bordering the narrow slip distribution holds a pointed light on the sequence evolution. The 2022 Adriatic compressional earthquake sequence confirms the ongoing seismotectonic activity of this sector of the Apennines that is still propagating toward its foreland roughly in a piggy-back thrust sequence.

Methods
We relocated 594 aftershocks in the period 22 November 2022-05 January 2023 and the 09 November mainshock by using the hypoDD relative method 8 . We used the P-and S-wave arrival times picked by the seismologists on duty in the INGV control room, obtaining 30,226 P-wave and 12,739 S-wave differential times computed from the initial P and S wave phase readings at 55 stations. We optimized the damping parameter and verified the convergence of the solution with a final rms of 0.13 s. The accuracy of the hypocenter solutions has been estimated by using a bootstrap approach. Starting from the residual vector derived from the double difference inversion, we use the bootstrap resampling method 36 to create 200 random residual vectors of unit weight that we add to the observed differential travel times. This results in 200 bootstrap datasets that are inverted as in the real case. For each event we thus obtain 200 values for the hypocenter parameters and this allows us to infer a realistic estimation of the uncertainties, that are within 100 m on average.
We report on the histograms of Fig. S5 of the supplementary material the average value -computed on 200 bootstrap samples-of the uncertainties on the X, Y, and Z coordinates. The mainshock relocation is obtained by using the hypoDD method after a re-picking of the P-and S-phase. In order to verify if that the hypoDD relocation does not cause a significant shift of the entire earthquakes cluster (a circumstance that may raise from the relative relocation method), we compare the hypoDD relocation shown in Fig. 3a, b with the 1D location. The result of this exercise is shown into a new Figure (Fig. S6) of the supplementary material.
To characterize the seismic sequence in terms of magnitude distribution and decay we have analyzed data from the seismic control room. By using as input the earthquake relocated catalog and by adopting a maximum likelihood estimation approach 37 , we estimated the magnitude distribution of the seismic sequence following the Gutenberg-Richter 23 recurrence relation and adopting: where a and b represent the annual level of seismicity and the ratio between the number of small and large earthquakes, respectively; N(M) is the cumulative number of earthquakes with magnitude M and larger within the used catalog.
Aftershocks/day and the daily average energy release is plotted as a function of time. The completeness magnitude threshold of 2.0 which is reached 4 h after the mainshock. Uncertainties on the rates were calculated by bootstrapping the data sets with 100 realizations and indicated with vertical bars at each data point in Fig. 3.
Raw daily GNSS observations were processed by using the GAMIT/GLOBK software 9,38 to estimate the coseismic displacements (additional details in supplementary material). The analyzed dataset includes 15 GNSS stations (Table 1) located within 60 km from the epicenter, from 1 October to 10 November 2022. Most of these stations are property of Eni S.p.A. and are located on onshore and offshore industrial infrastructures (e.g., seabedanchored platforms and storage centers); the other stations belong to the NetGEO network (http:// www. netgeo. it), owned by Topcon positioning and INGV RING network (http:// ring. gm. ingv. it ). To improve the overall configuration of the network and tie the regional measurements to an external global reference frame, data coming from 15 continuous stations belonging to EUREF (https:// epncb. oma. be), ASI (http:// geodaf. mt. asi. it) and FReDNet (https:// fredn et. crs. ogs. it/ DOI/) were introduced in the processing. The daily estimates of loosely constrained station coordinates have been combined in GLOBK to estimate a consistent set of daily coordinates www.nature.com/scientificreports/ (i.e., time-series) for all processed stations. Some stations show a significant offset on their time-series, therefore we computed the amount of 3D coseismic displacements, by differencing the average sites position in the three days before and the two days after the two main events with minimal constraints (i.e., constraining translations, scale, and rotations of the network solution to 0.1 mm). Due to the very short time between the occurrence of the two M w > 5 earthquakes, the estimated displacement pattern represents the cumulative effects of both shocks. Uncertainties associated with the coseismic displacements have been estimated as the "sum in quadrature" 39 of the errors associated with the average sites position before and after the events. As mentioned above, average uncertainty for the horizontal and vertical components are ~ 1 and 3.6 mm, respectively. These estimations must be considered as an overestimation of the real ones, due to the conservative approach used to their computation. Seismic source model is performed using the Okada formalism 26 , implemented in the SARscape software by SARmap (sarmap SA, 5.6 version), working in the ENVI® (L3HARRIS) environment. First, we performed an unconstrained non-linear inversion of the coseismic displacements measured by GNSS stations, to retrieve geometry and kinematics of the fault, followed by a linear inversion to estimate the slip distribution along the plane 40 . Because of the large uncertainties, the contribution of the vertical component in the uncertainty-weighted inversion is negligible, we considered the only horizontal vectors in the inversion. In supplementary materials we show the slip uncertainty (1σ predicted slip error standard deviation) along the fault plane (Fig. S7). To visualize the slip resolution in depth, we performed an additional linear model to retrieve the slip distribution along the fault plane using a variable size grid (Figs. S8 and S9). In our approach 41 patch dimensions are automatically determined by the inversion algorithm maximizing the model resolution matrix 42 . Patches vary from 3 × 4 km at shallow depths to 8 × 6 km at 15 km depth. Following the same approach we estimated uncertainties 40 .
The relationships between the mainshock source and the closest faults is investigated by using the Coulomb Failure Function. Results of seismic source modeling, seismological analysis and subsequent interpretation reveal the activation of a low angle thrust, corresponding to the ITCS106-Pesaro mare-Cornelia Composite Source, from DISS catalogue. Thus, we refined and tune the DISS source, adapting to our modelling result to calculate the stress increase along the activated fault plane. The underground fluid pressure is not considered and the friction coefficient is considered constant at 0.4. Being the absolute stress value unknown, the ΔCFF is defined as 43 : in which τ is the shear traction, σ is the normal traction (positive for traction) and μ' is the apparent friction coefficient. Stress change is computed assuming that the half-space is a Poisson solid (ν = 0.25) with Young modulus (E) of 75 GPa and a shear modulus (G) of 30 GPa.

Data availability
GNSS solutions are reported in Table 1  www.nature.com/scientificreports/ Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.